# Title     : TODO
# Objective : TODO
# Created by: Administrator
# Created on: 2019/7/24

library(optparse)
library(magrittr)
library(tidyverse)
library(ggrepel)
library(extrafont)
library(ggforestplot)
library(broom)

option_list <- list(
  make_option("--i", default = "AllMet_Raw.txt", type = "character", help = "metabolite data file"),
  make_option("--g", default = "group.txt", type = "character", help = "sample group file"),
  make_option("--cc", default = "calculate_config.json", type = "character", help = "config file"),
  make_option("--pc", default = "", type = "character", help = "config file")
)
opt <- parse_args(OptionParser(option_list = option_list))

df_nmr_results <- read_csv(
  file = "metabolite.csv"
) %>%
  mutate(Metabolite = factor(Metabolite, levels = unique(Metabolite))) %>%
  gather("SampleID", "Value", -Metabolite) %>%
  spread(Metabolite, "Value")

df_clinical_data <- read_csv(file = "clinical_data.csv") %>%
  rename(SampleID = Sample)

df_full_data <- left_join(
  x = df_nmr_results,
  y = df_clinical_data,
  by = "SampleID"
)

nmr_biomarkers <- df_nmr_results %>%
  select(-"SampleID") %>%
  colnames()

df_long <- df_full_data %>%
  select(all_of(nmr_biomarkers), gender,BMI) %>%
  mutate_at(
    .vars = vars(all_of(nmr_biomarkers)),
    .funs = ~.x %>% log1p() %>% scale %>% as.numeric()
  ) %>%
  gather(
    key = biomarkerid,
    value = biomarkervalue,
    all_of(nmr_biomarkers)
  )

df_long<-df_long %>%
  mutate(gender=factor(gender,levels=unique(gender)))

df_long

df_assoc_per_biomarker <- discovery_regression(
  df_long = df_long,
  model = "lm",
  formula =
    formula(
      biomarkervalue ~ BMI
    ),
  key = biomarkerid,
  predictor = BMI
)

df_assoc_per_biomarker

df_pca <- df_full_data %>%
  select(all_of(nmr_biomarkers)) %>%
  nest(data = everything()) %>%
  mutate(
    pca = map(data, ~prcomp(.x, center = TRUE, scale = TRUE)),
    pca_aug = map2(pca, data, ~broom::augment(.x, data = .y))
  )

df_pca_variance <- df_pca %>%
  unnest(pca_aug) %>%
  summarize_at(.vars = vars(starts_with(".fittedPC")), .funs = ~var(.x)) %>%
  gather(key = PC, value = variance) %>%
  mutate(cumvar = cumsum(variance / sum(variance)), PC = str_replace(PC, ".fitted", ""))

pc_99 <- df_pca_variance %>%
  filter(cumvar <= 0.99) %>%
  nrow()

psignif <- signif(0.05 / pc_99, 1)

df_to_plot <- df_assoc_per_biomarker %>%
  drop_na(.data$estimate)

df_to_plot <- df_to_plot %>%
  rename(name = biomarkerid)

df_to_plot <- df_to_plot %>%
  head(20)

forestplot(
  df = df_to_plot,
  pvalue = pvalue,
  psignif = psignif,
  xlab = "1-SD increment in biomarker concentration\nper unit increment in BMI",
)

nmr_biomarkers <- dplyr::intersect(
  ggforestplot::df_NG_biomarker_metadata$machine_readable_name,
  colnames(df_demo_metabolic_data)
)






